/*
  License...
*/

#include "morphology.h"

void Find(CImg* where_img, CPattern* what_img, double lowbound, CFindRes* res)
{
	double c0, c1, c2;
	double* projv = new double[what_img->size];
	ulong x, y, dx, dy, i, j, k;
    double funcval, funcnum, funcdenom, dd;
	res->numres = 0;

#ifdef SHOWPROGRESS
	ulong proc = (ulong) (where_img->width - what_img->width) / 20, procindex = 0, pp = 0;
#endif

	for(x = 0; x <= where_img->width - what_img->width; x++)
	{

#ifdef SHOWPROGRESS
		procindex++;
		if(procindex > proc)
		{
			procindex = 0;
			pp += 5;
			cout << pp << "% complete..." << endl;
		}
#endif

		for(y = 0; y <= where_img->height - what_img->height; y++)
		{
			c1 = 0;
			for(k = 0; k < what_img->len1; k++)
			{
				i = x + what_img->ind1[k].x;
				j = y + what_img->ind1[k].y;
				c1 += where_img->buff[ where_img->width * j + i ];
			}
			c1 /= what_img->len1;

			c2 = 0;
			for(k = 0; k < what_img->len2; k++)
			{
				i = x + what_img->ind2[k].x;
				j = y + what_img->ind2[k].y;
				c2 += where_img->buff[ where_img->width * j + i ];
			}
			c2 /= what_img->len2;

			for(k = 0; k < what_img->len1; k++)
			{
				i = what_img->ind1[k].x;
				j = what_img->ind1[k].y;
				projv[ what_img->width * j + i ] = c1;
			}

			for(k = 0; k < what_img->len2; k++)
			{
				i = what_img->ind2[k].x;
				j = what_img->ind2[k].y;
				projv[ what_img->width * j + i ] = c2;
			}

			c0 = 0;
			for(dx = 0; dx < what_img->width; dx++)
			{
				for(dy = 0; dy < what_img->height; dy++)
				{
					i = x + dx;
					j = y + dy;
					c0 += where_img->buff[j * where_img->width + i];
				}
			}
			c0 /= what_img->size;

            funcval = 0.0;
            funcnum = 0.0;
            funcdenom = 0.0;
			for(dx = 0; dx < what_img->width; dx++)
			{
				for(dy = 0; dy < what_img->height; dy++)
				{
					i = x + dx;
					j = y + dy;
					dd = where_img->buff[j * where_img->width + i] - projv[dy * what_img->width + dx];
					funcnum += dd * dd;
                    dd = projv[dy * what_img->width + dx] - c0;
					funcdenom += dd * dd;
				}
			}
			funcval = funcnum / funcdenom;

			if(funcval <= lowbound)
			{
				if(res->numres < res->buffsize)
				{
					res->res_dist[res->numres] = funcval;
					res->res_pos[res->numres].x = x;
					res->res_pos[res->numres].y = y;
					res->numres++;
				}
			}
		}
	}

	delete projv;
	return;
}

void Find2(CImg* where_img, CPattern* what_img, ulong num_where_img, double lowbound, CFindRes* res)
{
	for(ulong i = 0; i < num_where_img; i++)
	{
		Find(where_img, &what_img[i], lowbound, &res[i]);
	}
	return;
}